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Minimal models of self-propelled particles with short-range volume exclusion interactions have 
been shown to exhibit signatures of phase separation. Here I show that the observed interfacial sta¬ 
bility and fluctuations in motility-induced phase separations (MIPS) can be explained by modeling 
the microscopic dynamics of the active particles in the interfacial region. In addition, I demon¬ 
strate the validity of the Gibbs-Thomson relation in MIPS, which provides a functional relationship 
between the size of a condensed drop and its surrounding vapor concentration. As a result, the coars¬ 
ening dynamics of MIPS at vanishing supersaturation follows the classic Lifshitz-Slyozov scaling law 
at the late stage. 

PACS numbers: 


I. INTRODUCTION 


explain the scaling between the interfacial width and the 
system size recently observed in MIPS El- 


Phase separation is a ubiquitous phenomenon in nature 
and is manifested by the partitioning of the system into 
compartments with distinct properties, such as the differ¬ 
ent particle densities in the two co-existing phases in the 
case of liquid-vapour phase separation. Phase separation 
under equilibrium dynamics is a well investigated physi¬ 
cal phenomenon m- Recently, signatures of phase sep¬ 
aration have been reported in non-equilibrium systems 
consisting of active particles [Mg. It is therefore a nat¬ 
ural question to ask to what degree we can extend our 
knowledge of equilibrium phase separation to the phase 
separation phenomenon observed in active systems. In 
the case of minimal models of active particles with sim¬ 
ple volume exclusion interactions, the phenomenon of 
motility-induced phase separations (MIPS) has received 
considerable interest [3, 15H81-IT8]. In particular, the idea 
of an effective surface tension in motility-induced phase 
separations (MIPS) has been advocated US HE]. At the 
gas-liquid interface in thermal equilibrium, surface ten¬ 
sion results from the pulling of molecules at the interface 
due to their attractive interactions [20] , In a system of 
active particles with purely repulsive interactions, it is 
unclear how such “pulling” can occur as the particles 
can only push. To probe what happens at the interface, 
I study here the microscopic dynamics of the active par¬ 
ticles in the interfacial region by a combination of simu¬ 
lation and analytical methods. Specifically, using mean- 
field type arguments, I will demonstrate how pressure 
balance is achieved between the condensed phase and the 
dilute (vapor) phase, and how the Gibbs-Thomson rela¬ 
tion arises in a system where a circular condensed drop 
co-exists with the vapour phase. Furthermore, by incor¬ 
porating the stochastic nature of particle dynamics, I will 
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A. Motility-induced phase separation 


I will first focus on a minimal model system that ex¬ 
hibits MIPS in two dimensions (2D)— A collection of self- 
propelled particles with excluded area interactions that 
undergo rotational fluctuations. Specifically, the dynam¬ 
ical equations are 

^ Iu-IjD+^V, (1) 

= V2Dr9i(t) (2) 

where i is an integral index enumerating the particles 
in the system, Vj = cos 9iX + sin B^y with the angle Qi 
(with respect to the x-axis) being the orientation of the 
*-particle, gi(t) is a noise term with Gaussian probability 
distribution with zero mean and unit variance, D r sets 
the magnitude of the rotational fluctuations, U(.) cor¬ 
responds to the potential function for short-ranged area 
exclusion interactions, r) is the drag coefficient and f a is 
the constant active force that drives the particles in the 
system. In particular, u = fa/g is the constant speed 
of a particle when it is not within the area exclusion 
zone of another particle. Previous numerical work has 
indicated that phase separation in this minimal system 
occurs as u increases, but the actual form of U is unim¬ 
portant nun ini- For instance, U could be of the form 
of a Weeks-Chandler-Andersen potential [211 : 


U(r) 


A 

0 


(r) 12 — 2 (y) 6 ) + 1 


if r < a 
otherwise . 


(3) 


This will be the particular form of potential function em¬ 
ployed in this work. Also, the time and length units will 
be set by having a = 1 and D r = 3. Note that I will 
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focus exclusively on the non-equilibrium dynamics of the 
system and so translational Brownian motion is ignored. 


II. FLAT INTERFACE 


A. Point particles 


To understand the microscopic dynamics at the inter¬ 
face, it is instructive to first look at a system of point 
particles, i.e., the interaction potential U is zero. Even 
this simple system distinguishes itself from equilibrium 
system in that aggregation will spontaneously happen in 
the proximity of a force-absorbing but frictionless wall 
(left column of Fig. [I]). In other words, the particles are 
free to slide and rotate at the wall, but they cannot pen¬ 
etrate the wall. Further complex patterns are revealed 
when one looks at the particles’ orientation distribution 
as well as the position distribution (Fig. [lje)). At the 
wall, most of the particles are left going, as indicated by 
the high concentration of orientation at around 9 ~ tt. 
This results from the fact that only left-going particles 
remain at the wall. Just outside the wall, the distribution 
is highly peaked at 9 just below tt/2 and just above 37r/2, 
which reflects the particles’ orientation after they move 
away from the wall. The orientation anisotropy decays 
as one moves away from the wall. 

In this system, the pressure acting on the wall can be 
expressed as 


Pw = 



Xw{9) cos 9d9 


(4) 


where the orientation distribution function of the par¬ 
ticles at the wall per unit length is denoted by xwifi). 
Note that since we are dealing with a 2D system, the unit 
of pressure is [force] / [length]. 

To further analyse xw(0)i one can perform dimen¬ 
sional analysis to conclude that 

Xw{&) = , (5) 

V D r 

where P{9) is a function dependent only on 9 , and is 
the particle concentration far from the wall. To obtain 
the exact functional form oiJ-(O), one needs to solve a set 
of two coupled differential equations with mixed bound¬ 
ary conditions [32], whose solution consists of a series of 
Matliieu functions. Unfortunately, the expansion coeffi¬ 
cients in the series are not analytically tractable and so 
an analytical expression is lacking. However, P(9) can 
be readily estimated numerically (left column of Fig. [l]), 
which allows us to obtain the following: 


Pw 


f 2 o /' 37r/2 

/ T{9)cos9d9 
VPr J tt/2 


faPoo 

2r/D r 


( 6 ) 


The second expression is equivalent to the swim pressure 
of a system of active particles in 2D [35H3B] , which I have 
obtained numerically here. 


B. Repulsive particles 


Remarkably, much of what we have seen in the point 
particle case remains true when we add mutually repul¬ 
sive interactions to the particles. When a force absorb¬ 
ing but frictionless wall constitutes the left boundary of a 
semi-infinite system, phase separation occurs where the 
condensed phase is located close to the wall (Fig. [ijb)). 
Inside the condensed phase, the orientation is isotropic 
(Fig. [ljf )). The reason behind the isotropy is that the 
impeded motility of the particles render them staying put 
for a duration much longer than the orientation decoher¬ 
ence time ~ 1/D r . However, note that the particles’ 
locations are not frozen in time as shown by the black 
particles in Fig. [ljb) , which were the first column of 
particles next to the wall at the beginning of the sim¬ 
ulation. As one moves further to the right, one first 
encounter a layer of left-going particles (shown by the 
bright red patch centred at x ~ 17 in Fig. [ljf)). This 
represents the accumulation of particles with orientation 
highly centred at around 9 ~ 7r, analogous to the ac¬ 
cumulation of active point particles at the wall, except 
that the particles are now spread over a range of x posi¬ 
tions due to volume exclusion interactions. Further right¬ 
wards, we encounter the pattern of two escape trajecto¬ 
ries away from the interfacial region indicated by the two 
red-yellow branches emerging from the red patch. As one 
moves further away to the right, the orientation becomes 
isotropic again. From this discussion, it is clear that in 
the bulk of the condensed and vapour phases, the cor¬ 
responding orientation distributions are both isotropic, 
while in the interfacial region separating them there is a 
high level of orientation anisotropy. 

Let us now calculate of the force exerted on the wall by 
the active force of these particles. Since in the condensed 
phase, the orientation is isotropic, the pressure felt by 
the wall due to these active forces is 


o( a ) _ 

r w ~ 


&Pcfa 

j- 3tt/2 

/ cos 9d9 

2n 

Jtt/2 


a Pcfa 


(7) 


where p c is the concentration in the bulk condensed 
phase. Due to the orientation isotropy, the expression 
here scales like f a instead of f 2 in the point particles 
system (Eq. <§)• Besides the active force contribution, 
the wall will also feel additional forces arising from the 
repulsive interactions, which, we will see, constitute an 
important contribution in achieving pressure balance in 
the interfacial region in MIPS (Sect. IID). 


C. Locating the interface 

The concentration variation across the two phases 
shown in Fig. [2ja) is similar to typical equilibrium phase 
separation. What distinguishes MIPS is the high orien¬ 
tational anisotropy between the phases (Fig. §b) & (c)). 
As in equilibrium fluids, the location of a sharp interface 
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FIG. 1: Steady-state configurations of active particles confined by a force absorbing wall on the left: point particles (left 
column) and repulsive particles (right column), a) & b): A snapshot of the system at the end of the simulations with the 
orientations depicted by the blue arrows. The red circles depicted in (b) are of diameter a = 1. The wall at x = 0 is perfectly 
force absorbing (see A 1 for simulation details), c) & d): The histograms show the horizontal distributions of the particles, e 
& f): The colourmaps show the deviation from the mean in the particles’ orientations at different horizontal positions. The 
colour scale corresponds to the measure: 2 (hi(x) — (hi(x)))/ max^ hi(x) where i is the row index and hi(x) is the frequency. 
The simulation parameters are: f a = 100, y = 1, a = 1, D r = 3, A = 25/6. 


between the two phases can be defined somewhat arbi¬ 
trarily m ■ In our case, since the pronounced minimum 
of ( v x ) = (cos 6) is easy to locate (indicated by the red 
broken line) and its location also marks the onset of the 
increase in Q yy = — (cos(2#)/2) (the yy component of 
the nematic order parameter Q) [2S], which signifies the 
escape of particles from the condensed phase, it is a con¬ 
venient choice for the interface location. In other words, 
this convention implies that right outside the interface of 
the condensed phase, there is a layer of particles trav¬ 
elling preferentially along the interface, as indicated by 
the peak in Q yy (Fig. [2](c)) . The active forces of these es¬ 
caped particles are potentially the cause of the emergence 
of a negative surface tension according to its mechanical 
definition [18). The definition of the interface location 
has of course no physical significance, but this does pro¬ 
vide a working definition useful for the sharp interface 
model discussed below. 


D. Interface stability 

1. Pressure balance at a sharp interface. 

In this section we will see how pressure balance can be 
achieved in MIPS. Note that the discussion in this sec¬ 
tion amounts to a simplified exposition of that in |23ll29| . 
Its presentation here is for the self-containedness of the 
paper and will help us understand the approximations 
used in later sections 

I will start by discussing a sharp interface model. In 
this drastically simplified model, let us imagine that the 
phase separated system is partitioned by a sharp inter¬ 
face where in the vapour phase, the concentration is low 
enough that the system behaves like a system of active 
point particles, and in the condensed phase, the orien¬ 
tation distribution is isotropic. One can imagine such a 
system by first rotating Fig. [lje) by 180° and then col¬ 
lating it to Fig. [ljf) on the left. As calculated before, 
The pressure exerted by the vapour phase on the sharp 
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X 


FIG. 2: a) Particle concentration as a function of x. Same fig¬ 
ure as in Fig. 0d)- b) The mean horizontal component of the 
particles’ orientations ( v x ) vs. x (+ symbols, left y-axis), and 
pd x p vs. x (o symbols, right axis), c) The yy component of the 
nematic order parameter Q vs. x, where Q yy = — (cos(20)/2) 
so that high Q yy signifies that the orientations of the particles 
are preferentially pointing up or down. Similar to equilibrium 
fluids, defining the location of a sharp interface is somewhat 
arbitrary [271. The working definition proposed here is that 
the interface is set to be at the pronounced minimum of | (v x ) \. 

interface is f^p v /(2riD r ) (Eq. ©), this pressure is bal¬ 
anced by the pressure exerted by the condensed phase: 
apcfa/n + Pr, where the first term comes from the active 
force (Eq. 0 ) and the second denotes the pressure arisen 
from the repulsive force due to the area exclusion inter¬ 
actions. In other words, pressure balance is achieved 
if 


quired for MIPS. I assume for simplicity that at the onset, 
i) the condensed phase is not very compressed and so we 
can ignore P r , m and ii) the concentration ratio between 
the condensed phase and the dilute phase (p c /p v ) is of 
order 1, which let us to the minimal force requirement 
below for MIPS: 

fmin > ZapPr . ^ 

7T 

The above condition comes from the pressure balance 
condition at the interface alone. Interestingly, Eq. ([9]) 
reproduces the same scaling as obtained by Redner et al 
via a different approximation m- The result also sup¬ 
ports the notion that the Peclet number (Pe), usually 
defined as Pe oc f a /(ariD r ) in this context, is a key con¬ 
trol parameter in MIPS [TT1 [5CTj . 

From this sharp interface model, we can now see why 
the condensed phase with a high density of active and re¬ 
pulsive particles can remain stable against a backdrop of 
dilute concentration of active particles - the active parti¬ 
cles in the vapour phase impact an active pressure that 
scales as /„ directed towards the normal of the interface, 
while the countering active pressure from the condensed 
phase scales as f a ■ The quadratic dependence in the 
active force comes from the fact only particles pushing 
against the interface will remain on the interface while 
particles with orientation away from the interface will 
leave. The escapes of these particles thus open up space 
for yet other particles that serve to push against the in¬ 
terface. Form this perspective, the low concentration in 
the vapour phase is paramount for the stability of MIPS, 
for otherwise the particles with orientations away from 
the interface may be unable to leave effectively. 


2. Force balance in an interface of finite width. 


aPcfg p _ Pvfg 

7r r 2pD r ' 


( 8 ) 


For the simulation parameters used in Fig. |T] (with the 
units set by a = 1 and D r = 3), ap c f a /ir ~ 30, P r ~ 230 
and p v f%/(2r)D r ) ~ 330. We thus see that in the simu¬ 
lated system, the L.H.S. and the R.H.S. of Eq. ([s]) are of 
the same order of magnitude, indicating that the pressure 
balance condition is qualitatively satisfied. In addition, 
we see that much of the active force coming from the 
vapour phase is used to compress the condensed phase 
via the pressure P r . Note that although Eq. Q provides 
a pressure balance condition for the system, it does not 
mean that any system satisfying this condition is stable 
as it may still be unstable against fluctuations. This is 
not dissimilar to equilibrium fluids where pressure bal¬ 
ance together with chemical potential balance are needed 
to achieve phase stability. 

The pressure balance condition in Eq. ([8]) already al¬ 
lows us to estimate crudely the minimal active force re¬ 


Let us now go beyond the previous sharp interface pic¬ 
ture and see what happens within the interfacial region 
from the view point of particle dynamics. Ignoring fluc¬ 
tuations, the stability of the interface means that if a 
particle happens to be lying at the interface will, on av¬ 
erage, remain put. In our minimal model, a particle can 
only move due to two reason: i) its own active force driv¬ 
ing it to move in the direction dictated by its orientation, 
ii) repulsive force that pushes it away from its neighbours 
if it is of less than unit distance away from them. While 
the second force is common in both active and passive 
(equilibrium) systems, the active force is unique to non¬ 
equilibrium systems. Consider now a particle located at 
xq inside the interfacial region, i.e., where (v x (xo)} is 
varying (Fig. [2j. Since {v x (x 0 )) < 0, the active force 
will on average drive this particle to the left. On top of 
this, there are repulsive forces coming from neighbouring 
particles on the right hand side f r (X q + Ax). For the 
particle to remain still, the sum of these forces has to be 
countered by the repulsive forces coming from the left. 
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Therefore, 


A. Gibbs-Thomson relation: dimensional analysis 


fr{x 0 - a/2) = fa(v x {x o)) + f r (x o + a/2 ) . ( 10 ) 

Since the repulsive forces come from the repulsive poten¬ 
tial function U, let us replace the repulsive force by the 
pressure P p (x ) (the subscript p for passive) resulting from 
the corresponding system with the same particle config¬ 
uration and interaction potentials, but with the active 
force omitted. Since f r {x) — aP p (x), Eq. (10) leads to 


fa{v x {x o)) ^ a[P p (x o + a/2) - P p (x 0 - a/2)] (11) 

= . ( 12 ) 

da; 

In principle, P p {x) depends on the exact configuration of 
particles in the system, but if one adopts the simplifying 
assumption that the passive pressure depends solely on 
the particle concentration, one can then expand P p (x) 
with respect to the concentration p(x): 

P p {x) = c 0 + C!p(x) + c 2 p(x) 2 + 0(p 3 ) . (13) 


For equilibrium fluids, this is of course the virial expan¬ 
sion where Cq = 0 and C\ = k R T [27]. Since our system is 
fundamentally non-equilibrium (no translational Brown¬ 
ian motion, i.e., /cgT = 0 ), there is no guarantee that the 
same would apply here. But let us assume that such an 
expansion is possible in our system, then since P p comes 
purely from the repulsive interactions between particles. 
We thus expect that cq = 0 = c\ because as the con¬ 
centration goes to zero, there would not be any pairwise 
interactions. Therefore, the first non-trivial term in the 
expansion is c 2 p 2 . Note that c 2 > 0 since P p arises 
purely from the repulsive interactions. One could also 
incorporate active pressure into the analysis as, for ex¬ 
ample, done by Winkler et al 


Here, to order 0(p 2 ), Eq. (12) then leads to 


fa(v x {x o)) = 2a 2 c 2 p(x) 


dp(zo) 

da; 


(14) 


Remarkably, simulation result shown in Fig. [2jb) indeed 
seems to vindicate Eq. (14). 


III. CIRCULAR INTERFACE 


We have seen in the previous section how pressure bal¬ 
ance is achieved at a flat interface. However, previous 2D 
simulation studies have shown that similar to equilibrium 
phase separation, if the condensed phase in MIPS does 
not span the system size, the condensed phase is circu¬ 
lar. Here, we will see how the curvature of the interface 
affects the particle dynamics at the interface, and its con¬ 
sequence in terms of the coarsening dynamics. We will 
first study the emergence of the Gibbs-Thomson relation 
by dimensional analysis. 


In equilibrium phase separation, the Gibbs-Thomson 
(GT) relation dictates that the concentration (f> R right 
outside a droplet (of the condensed phase) of radius R is 

0R = (f + , (15) 

where </>q is the supersaturation concentration, i.e., the 
threshold concentration beyond which phase separation 
occurs, and v = is the capillary length with 7 be¬ 
ing the surface tension and v being the volume of the 
molecule. Since the concentration in the vapour phase 
outside a big drop is lower than that outside a small 
drop, a diffusive flux is set up that transfers material 
from the small droplet to the big droplet. This is the 
Ostwald ripening mechanism that dominates the phase 
separation kinetics at the late stage for systems with a 
small supersaturation [31]. 

I will now discuss why the GT relation would arise nat¬ 
urally in our active system. In the minimal MIPS system 
considered, the only parameters in the dynamical equa¬ 
tions are the free roaming speed u = f / 77 , the rotational 
diffusion coefficient D r: and the length scale of the short 
range area exclusion interaction a. Denoting now the 
vapour density far from a drop of radius R by p* R) and 
the density inside the drop (the condensed phase) by p c , 
then by dimensional analysis we have 


p c \D r a D r R 


(16) 


where T is some unknown scaling function dependent 
on its two dimensionless arguments. If we now assume 
that J- is regular with respect to the second argument in 
the sense that a Taylor series expansion exists (around 
u/D r R = 0), then the ratio above can be re-expressed as 


* 

Pr 

Pc 


H + K 


u 

D r R 


+ 0 



(17) 


where H and K are now just dimensionless functions 
of the first argument ( u/D r a ), i.e., R- independent. In 
terms of Eq. can be re-written as 



where v = -jf/y -, which may be termed the effective cap¬ 
illary length. In the large R limit, Eq. (18) becomes 
exactly the GT relation in Eq. (15). This analysis pro¬ 
vides an intuitive reason why one would naturally expect 
the GT relation to emerge as the drop radius grows in 
MIPS. 


B. Gibbs-Thomson relation: numerics 


I will now test Eq. (18) by simulating a coarse-grained 


model of MIPS. Let us first consider what happens to 
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FIG. 3: a) A schematic of the interfacial condition at a curved 
interface of curvature R~ . The particle is assumed to occupy 
a zone of diameter a and the particle can leave the droplet if 
its orientation is within the escape range of 2o\r indicated, b) 
A schematic showing a droplet (blue) in the condensed phase 
of radius R located at the origin co-existing with the dilute 
medium (vapour phase). An active particle (red circle) with 
orientation 9 (blue arrow) is located at the position (r, <p). 
The angle if> equals the difference between the orientation 9 
and the azimuthal coordinate ip. c) A snapshot of a simulated 
system with 1000 active point particles (red dots with orien¬ 
tations indicated by blue arrows) in an annular system with 
inner radius R = 100 and outer radius R + L r = 150. 



FIG. 4: a) The variation of the vapour concentration Pr{t) = 
- f dt/j(R(r,ip) away from the interface of a droplet of sizes 
R = 20, 60,100. The concentration decays rapidly from the 
interface and reaches a stationary value p* R a short distance 
away from the interface, b) Pr/pIo vs. R. The curve decays 
to 1 like R~ 1 , which is consistent with the Gibbs-Thomson 
relation (Eq. (181). Blues crosses are the simulation results 
and the red curve depicts the function 1 + 1.84iZ _1 . Note 
that the constant p* x is estimated numerically from the Gibbs- 
Thomson relation. See IA 21 for simulation details. 


an active particle in the condensed phase at the periph¬ 
ery that is curved (Fig. [3j)a)) . For such a particle, I as¬ 
sume that it occupies a zone of radius a (shown in red) 
sandwiched by two zones occupied by two neighbouring 
particles (light blue). Note that since the concentration 
at the interface may not reach the level of optimal pack¬ 
ing concentration (~ 0.91), a should be greater than the 
particle’s diameter a. Indeed, Fig. [2ja) suggests that 
the concentration is around 0.72 at the interface, which 
indicates that a ~ 1.2. To incorporate the effects of the 
neighbouring particles on the pink particle, I assume that 
as a result of the caging effect, the particle can only move 
out of the droplet if its orientation is within the 2er range 
depicted. Based on the diagram shown in Fig. [3j a simple 
trigonometric exercise leads to 

+ H + ■ (w) 

As expected, a lower curvature leads to a smaller escape 
range (smaller ctr). 

To analyse how the variation in the escape orienta¬ 
tion range affects the phase separated system at the 
steady-state, I consider a system with one condensed 
drop of radius R co-existing with the vapour phase (Fig. 
§b))- Let me denote the particle distribution func¬ 
tion in the vapour phase by pn(r,tp,d) where the first 
two arguments correspond to the particles’ locations and 
the last argument to the particles’ orientations. Due 
to rotational symmetry, one can eliminate one angu¬ 
lar argument by introducing the variable ip = 9 — tp 
[35] , and study instead the reduced distribution function 
(pn(r,ip) = f 0 pn{r, ip, (p + ip) r cos tpdip. On the drop’s 
periphery, the corresponding reduced distribution func¬ 
tion is denoted by Xr(V 0- Since the periphery is assumed 
to be infinitely thin, \R is only a function of ip and hence 
dimensionless. In addition, I assume that drops of all 
sizes have the same interior concentration p c , xr(V ; ) is 
thus related to the p c as follows: J dV’Xfl(V’) — 2 dRp c . 

To study the distribution functions, I assume again 
that the vapour concentration is low enough that pairwise 
repulsive interactions can be ignored, and simulate the 
dynamics of active point particles, i.e., non-interacting 
active particles, in an annular geometry of inner radius 
R and outer radius R + L r ( Fig. [3](b)). As in the linear 
case, the particles’ orientations are randomised if they 
reach the outer circular boundary, while if the i-tli parti¬ 
cle reaches the inner boundary, its positions will remain 
fixed until its orientation is within the escape range, i.e., 
until ipi is between — or and ctr. Simulation results are 
shown in (Fig. [4]). Away from the interface, it is observed 
that the concentration rapidly reaches a stationary value 
P*r=\ f d i’CRir, ip) for, say, r > R + 20 (Fig. 0a)). As 
expected from previous discussion, the vapor concentra¬ 
tion p* R goes down with R since a flatter interface leads 
to a narrower escape range, which leads to a smaller out- 
flux of particles from the condensed phase. Fig. ®b) 
shows that p* R decays to p* x (the vapour concentration 
as R — > oo) like i? -1 , which, as we have seen, is consistent 
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FIG. 5: A schematic depicting a wavy interface where the con¬ 
densed phase is depicted in blue. The location of the interface 
(purple) is given by the function h(y). Due to the caging effect 
from neighbouring particles, the red particle at the interface 
will have a higher chance of escaping compared to the green 
particle because the escape orientation range (grey area) is 
bigger. 


the Gibbs-Thomson relation in equilibrium systems. 


IV. FLUCTUATING INTERFACE 


I have so far ignored fluctuations in the interfacial pro¬ 
file. In reality, the interface of course fluctuates, which is 
already discernible from the spatially constrained system 
shown in Fig. 0 b). In particular, previous simulation 
result points to the scaling law [IS] : 

w 2 l ~L, (20) 

where wl is the steady state interfacial width: 

WL y~J^\jj' 0 (%) 2 - h 2 )dy (21) 


with h being the average position of the interface. Here, 
the symbol h(y) denotes the location of the interface, 
i.e., the location of the peak of (v x ) (Fig. [2]). To under¬ 
stand the scaling observed, let us consider the effects of 
interface curvature on the particle exchange dynamics. 
Although the previous section focuses only on a circular 
interface, i.e., the curvature is positive, one can easily 
extends Eq. (19) to allow for concave interface as well 
(Fig- 0- The physical motivation behind the formula is 
the same, a particle at a highly convex portion of the 
interface will have a wider escape orientation range (red 
particle in Fig. [5]) than a particle at a highly concave 
interface (green particle). 

Since the fluctuations ultimately come from the fluctu¬ 
ating dynamics of particle exchange at the interface, we 
need to model the steady-state dynamics of h, stochasti¬ 
cally. The simplest equation of motion (EOM) for the 
interface that incorporates both the effects of stochastic- 



FIG. 6 : Interface fluctuations as measured by the interface 
width wLit) = X)f=i (hi{t) — h(t)) 2 . The curve collapse of 
systems with difference linear dimension L upon rescaling is as 
predicted by the EW model. Inset plot: The interface width 
at the final time iut(t/) shows a linear dependence with \fh 
where L is the system size. Blues crosses are the simulation 
results and the red line is a guide for the eyes. See Appendix 
I A 31 for simulation details. 


ity and curvature-modified outflux is 


— = 2a6 1 (y) - 2a'(1 + j3n{y))b 2 {y ) , (22) 

where a denotes the rate of particles coming into the in¬ 
terface, and thus contributing to the growth of h, while 
a'(l + (i'K{y)) denotes the rate of particle escaping with 
the effect of local interface curvature (n(y) = d 2 h/dy 2 ) 
taken into account. The noise terms are bi(x) which are 
Markovian, spatially independent and are either 0 or 1 
with equal probability. Since the interface does not move 
at the steady state by assumption, a has to be identical 
to a'. From now on, I will focus exclusively on the hydro- 
dynamic limits (large temporal and spatial scales). So let 
us coarse grain h by defining a new coarse-grained height 
function h(y ): 


1 rv+t /2 _ 

%) = 7 / %')<V (23) 

t Jy-l /2 

where a -C i -C L y and I is large enough that 
Iy-i /2 bi(y')<ly' become Gaussian as a result of the cen¬ 
tral limit theorem. The EOM of h(y) is then 


dh 

dt 


a 


a 


9i + 2 ~ (1 + P K ) (^92 + 7 ^ 

f3 K 

— + Png2 + (<7i - 92) 


(24) 

(25) 


where gi are now Gaussian noises such that 


(9i(y,t)) = 0 
( 9 i(y,t)gj(y’,t ')} = 


(26) 
y')- (27) 



















In the long-wavelength limit, the fluctuating term 
a(3ng2 ~ d 2 h/dy 2 —» 0 and so the only relevant fluc¬ 
tuations come from the Guassian fluctuations a(gi — < 72 )- 
Therefore, in the hydrodynamic limits, Eq. ( |25[ ) is ex¬ 
actly the Edwards-Wilkinson model [Mj, with the effec¬ 
tive surface tension given by a/3/2. Note that the ef¬ 
fective surface tension here is always positive, which is a 
requirement for having a stable interface. As such this 
effective surface tension is distinct from the mechanical 
definition of the surface tension, which has been shown 
to be negative in MIPS [IS]. As mentioned in Section 
2.2, the negative tension around the interface is likely to 
be caused by the particles escaping from the condensed 
phase that are now travelling close to being parallel to 
the interface. Consistent with our definition of the loca¬ 
tion of the interface (Fig. [2]), these escaped particles are 
not considered part of the condensed phase and are thus 
ignored in our discussion of the interfacial fluctuations. 
In other words, the effective surface tension derived here 
is distinct from the mechanical definition of the surface 
tension discussed by Bialke et al [18] . 

Given that our stochastic model is equivalent to the 
Edwards-Wilkinson model, the temporal and steady- 
state dynamics of interfacial width is known to following 
the scaling form: 


w L (t) = L a J ew ( 28 ) 

for some scaling function J r EW (-) (Fig. [6). In a 2D system 
where the interface is a line, a = 1/2 and /? = 1/4 [MI PT5l. 
As shown in Fig. [6] these expectations are confirmed with 
direct simulation of a discretised version of the original 
EOM in Eq. ( [22] ) . This model thus provides an analytical 
argument supporting the steady state scaling wi,{t —► 
00 ) ~ L 1 / 2 recently observed numerically [TS] . 

To summarise this section, I have incorporated the 
caging effect as discussed in Sect. m into the modelling 
of the stochastic dynamics of particle exchanges at the 
interface. The model equation is then shown to be equiv¬ 
alent to the Edward-Wilkinson model in the hydrody¬ 
namic limits. In particular, the emergence of the effective 
surface tension term (a/3n/2) from the particle dynam¬ 
ics at the interface also explains why the interface is flat 
when both phases span the system, and circular when 
one phase does not span the system. 


V. MIPS IN 3D 


I have so far analysed the interfacial properties in 
MIPS in 2D using a combination of analytical and nu¬ 
merical methods. Here I will extrapolate the results ob¬ 
tained to MIPS in 3D. 


A. Flat interface 


Employing the sharp interface model for MIPS in 3D, 
the pressure balance equation in Eq. (|8| becomes 


aPcfg _ Pv fg 
2 tt ’ 6 yD r 


(29) 


where on the R.H.S. the active force contribution cor¬ 
responds to the swim pressure of active particles in 3D 
in the vapour phase [23J, and on the L.H.S., the first 
terms comes from active contribution to the force as¬ 
suming again that the orientation is isotropic in the con¬ 
densed phase. With regards to the minimal active force 
required for MIPS, using the approximations that f r is 
negligible and that p c /p v — 1, we arrive at 


rmin 
J a — 


3arjD r 


(30) 


which is very similar to the expression in 2D (Eq. ©)• 


B. Spherical interface 


The dimensional analysis presented in Sect. Ill A ap¬ 
plies also spherical drops in 3D. Therefore, if the assump¬ 
tion that the escape range decreases with the curvature 
of the drop, then the Gibbs-Thomson relation should 
emerge in the large drop limit (Eq. (]l8|). In particu¬ 
lar, we again expect the MIPS coarsening kinetics to be 
equivalent to the equilibrium scheme at low supersatura¬ 
tion [1511221125]- 


C. Interface fluctuations 

For MIPS in 3D, the interface is two dimensional and 
so there are two principal curvatures. If we adopt the 
natural assumption that the escape range now depends 
on the mean curvature, then the theoretical analysis pre¬ 
sented in Sect. |IV| applies in 3D straightforwardly. As a 
result, keeping only the linear terms will again lead to the 
Edwards-Wilkinson model in the hydrodynamic limits. 


VI. DISCUSSION & OUTLOOK 

In this paper I have investigated the microscopic dy¬ 
namics of active particles in the interfacial regions in 
MIPS using a combination of simulations and analytical 
arguments, and demonstrated i) how interface stability is 
achieved, ii) why the GT relation emerges in MIPS, and 
iii) how interface fluctuations scale with the system size. 
Therefore, I have shown that all the observed “surface 
tension” related phenomena found in MIPS result from 
the microscopic dynamics of the active particles. More 
specifically, I have demonstrated that pressure balance in 
MIPS is achieved because of the orientation anisotropy in 
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the region, which leads to a high active force directed to¬ 
wards the condensed phase. By incorporating the caging 
effects of neighbouring particles in the peripheral of a 
condensed drop, I have shown how the Gibbs-Thomson 
relation emerges naturally in MIPS, which dictates that 
the larger the condensed drop, the smaller the vapour 
concentration outside the drop. If the supersaturation 
level is small, the GT relation leads to diffusive transfer 
of active particles from small drops to larger drops. As 
a result, the late-stage coarsening kinetics in an active 
phase-separating system should follow the temporal scal¬ 
ing as in equilibrium phase separation: i.e., the average 
droplet size in the system (. R(t )) goes like f 1 / 3 [To 324201 . 
In addition, the droplet size distribution should approach 
asymptotically the universal size distribution obtained 
by Lifshitz and Slyozov [32]. Lastly, motivated by the 
same caging effects, I have proposed a stochastic model 
that describes the interfacial fluctuations in MIPS. In this 
model, the probability of particles leaving the interface 
is assumed to be dependent on the interface curvature. 
Analytical argument is then provided to show that the 
proposed model belongs to the same universality class of 
the Edwards-Wilkinson model. 

There are a number of future directions that are of in¬ 
terest. For instance, phase separation may play a role in 
the cytoplasmic re-organisation during asymmetric cell 
division [371 How the activity in the cytoplasm 

due to the many motor proteins contribute to such re¬ 
organisation via phase separation awaits more attention. 
Moreover, the fact that active phase separation occurs 
naturally begs the question of the existence of the critical 
point as in the equilibrium case. The critical transition 
in incompressible active fluids has recently been shown 
to give rise to a novel universality class [35]. And if a 
critical transition exists in MIPS, will the critical expo¬ 
nents be identical to those in the equilibrium case which 
belong to the 2D Ising Universality class? This question 
awaits further investigation. 


Appendix 


Appendix A: Simulation details 


boundary condition is enforced in the y, direction; while 
in the x direction, if x* +At < 0, then it is reset to zero, 
and if x* +At > L x , then x\ +At = L x and f?- +Ai is an 
angle chosen at random. 

The parameters of the simulations are: At = 10 -5 , 
a = 1, r) = 1, f a = 100, D r = 3, A = 25/6 for repul¬ 
sive particles and A = 0 for point particles. The system 
has width L x = 50 and height L y = 10sin(7r/3), with 
300 particles initialized in the configuration of a hexago¬ 
nal lattice (with spacing 1) next to the left boundary, and 
random orientations. Two hundred million time steps are 
evolved to equilibrate the system and then data are col¬ 
lected in the subsequent two hundred million time steps. 


2. Point particles in an annular geometry 


For the simulation results reported in Sect. Ill (Figs[3j) , 
the system is now annular with inner radius R and outer 
radius R + L r . The same updates as in Eqs (All and 
(A2) are performed for the point particles in the bulk 
of the system. Concerning the boundary conditions, if 
r l < R then the particle becomes part of the condensed 
drop periphery, and the orientation follows the update: 


g t+At = e t + ^2D^tgl , (A3) 

while the position remains the same until ip* = |0 f — 
¥>*1 < °R ( F ig-[lb)), in which case the particle leaves the 
condensed phase. Ifr‘ +At > R+L r , then r - +At = R+L r 
and d\ +At is an angle chosen at random. 

For distinct annulus geometry, the density pr{t) is nor¬ 
malised by p c , which is the density of particles at the 
inner circular wall. This is based on the assumption that 
p c is the same irrespective of drop sizes. 

The parameters of the simulations are: At = 10 —3 , 
t] = 1, f a = 100, D r = 3, A = 0. The system has 
1000 particles initialized with random orientations and 
positions. Twenty million time steps are evolved to equi¬ 
librate the system and then data are collected in the sub¬ 
sequent twenty million time steps. Simulations are per¬ 
formed for R = 20,40,60, 80 and 100 while L r is always 
50. 


1. Particle dynamics simulation 

For the simulation results reported in Sect. HU (Figs [I] 
& [2| , I numerically integrating the Langevin equations 
for each particle in the bulk of the system of size L x x L y 
by using the following updates: 

0t +At = ei + y/2^Atg\ (Al) 

= -))E V r^(l r ‘- r ‘l) + V C ° S ^ ’ (A2) 

where g\ are Gaussian distributed random variables with 
zero mean and unit variance, U is given by the Weeks- 
Chandler-Andersen potential shown in Eq. (|3|. Periodic 


3. Fluctuating interfaces 


To simulate interface fluctuations according to Eq. 


(22), I discretise the interface (a line) into L sites with 
height values hi where i = 1,..., L. Periodic boundary 
condition is enforced. The updates are performed as fol¬ 
lows: 


h^ At = hl + [l + P(h t i _ 1 + h\ +l — 2h t i )\g\ 

-[1 - P{h\_ 1 + h\ +1 - 2h\)]g t i , (A4) 

where g\ are either 0 or 1 chosen with equal proba¬ 
bility, and /3 = 0.1. The system sizes simulated are 
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L = 160,200,240,280,320 and 360. The total number 
of time steps simulated for each system size is 5L 2 /2. 
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